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Abstract 

We present a mesoscopic model for reactive shock waves, which extends the model proposed 
in [9]. A complex molecule (or a group of molecules) is replaced by a single mesoparticle, 
evolving according to some Dissipative Particle Dynamics. Chemical reactions can be handled 
in a mean way by considering an additional variable per particle describing the progress of 
the reaction. The evolution of the progress variable is governed by the kinetics of a reversible 
exothermic reaction. Numerical results give profiles in qualitative agreement with all-atom 
studies. 

1 Introduction 

Whereas multimillion atoms simulations are nowadays common in molecular dynamics studies 
with simple potentials, the time and space scales numerically tractable are still far from being 
macroscopic. The situation is even worse for nonequilibrium processes such as shock and 
detonation waves. Indeed, the simulation of detonation requires the description of a thin shock 
front, moving at high velocity, usually using a complicated empirical potential able to treat 
chemical events happening (dissociation, recombination - see [T] for a fundamental reference). 
To this end, toy molecular models were proposed at the early stages of molecular simulations 
of detonation (see e.g. [2]), until the first all-atom studies in the nineties [3] 2], allowed 
by the development of bond order potentials. Nevertheless, despite recent advances in the 
development of reactive force fields [5], the simulation of the decomposition of real materials 
still remains confined to small systems sizes and short time scales (see for example [6] for a 
state of the art study). A direct simulation of the macroscopic reaction zone of a real material 
is well beyond current computer capabilities, highlighting the importance of the development 
of reduced model for detonation. 

Following the pioneering works of Holian and Strachan [8], a reduced model for (inert) 
shock waves was proposed in [9], where complex molecules are replaced by mesoparticles. 
These mesoparticles are described by their positions, momenta, and have an additional degree 
of freedom: their internal energies. This model is strongly inspired by Dissipative Particle 
Dynamics (DPD) [10] with conserved energy [HIE]. In reduced models for shock waves [HIE], 
one mesoparticle stands for one complex molecule. Reduced models are interesting in practice 
to simulate large systems, and as an intermediate step in a truly multiscale approach, where 
some parts of the system could be treated with all-atom models, while the remaining parts 
would be treated with hydrodynamic models, implemented with particle discretizations such 
as Smoothed Particle Hydrodynamics |13l I14j . A first step to such a general formalism in the 
equilibrium case is proposed in [15] . 



In the reactive case, exothermic chemical reactions are triggered when the shock passes, 
and the energy liberated sustains the shock. To model detonation at the mesoscopic level, we 
introduce a new variable per mesoparticle, namely a progress variable, which characterizes the 
progress of the chemical decomposition. The dynamics can then be split into three elementary 
physical processes: 

(i) the translational dynamics of the particles, given by the dynamics of inert materials [9]; 

(ii) the evolution of chemical reactions through some kinetics on the progress variable; 

(iii) the exothermicity of the reaction: energy transfers between chemical and mechanical 
plus internal energies have to be precised. 

The paper is organized as follows. After summarizing the dynamics for inert materials, 
we turn to the evolution of the progress variable and the treatment of the exothermicity 
in the reactive case. A numerical implementation relying on some splitting based on the 
decomposition of the dynamics into elementary physical processes is also proposed. Numerical 
results eventually confirm the correct behavior of the model. Let us emphasize that we aim 
here at proposing a dynamics in qualitative agreement with all-atom detonation simulations, 
giving correct orders of magnitude for the speed of the detonation front and the width of the 
reaction zone. This dynamics is of course parametrized by a certain collection of parameters 
describing the initial and the final material, as well as the chemical kinetics. We however tried 
to limit their total number to keep only the essential ones, and give some tracks to estimate 
these parameters from small all-atom numerical simulations or physical experiments. 



2 A reduced model for inert shock waves 

We briefly recall here the inert model presented in [9]. We consider a system of N meso- 
scopic particles with positions (qi, .. . ,gjv), momenta (pi, . . . ,pjv), masses (mi, . . . ,mjv), and 
internal energies (ei, . . . , ejv). The internal temperature Tj is defined as T^ 1 = ds j^' ' . The 
function s = s(e) is the microscopic equation of state of the system, and relates the micro- 
scopic entropy a (arising from the existence of degrees of freedom not explicitely represented) 
and the internal energy e. As a first approximation, s(e) = C v \ne. Denoting by T the refer- 
ence temperature, (3 = l/(knT), and V the interaction potential, the system evolves in the 
inert case as 

dqi = — dt, 
rrii 

dpi = ^ -V gi Vjrij) dt - jjjX 2 (nj)Vij dt 

+*x(r ij )dW ij , (1) 

where r%j = \q% — qj |, Vij = — x is a weighting function (with support in [0, r c ] for some 
cut-off radius r c ), and the processes Wij are independent d-dimensionnal standard Wiener 
processes such that Wij = —Wji, The fluctuation-dissipation relation relating the frictions 
7ij and the magnitude a of the random forcing is 

Notice that we implicitely defined the temperature Tij as the average internal temperature: 
Tij = (Ti + Tj)/2. It can be shown that the measure 

dp(q,p, e) = Z~ 1 e~ l3H( - q ' p) exp -jj^ - /3ej \ 5 E =E 8p=p dqdpde 

is an invariant measure of the dynamics. In practice, we set a 2 = 2y/f3 and therefore 7ij/7 = 
T /Tij. The parameters used in the inert dynamics can be estimated as in [9]. 
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3 Evolution of the progress variable 



In the reactive case, chemical reactions are triggered when the shock passes. To model the 
progress of the reaction, an additional degree of freedom, a progress variable Ai, is attached 
to each particle. For the model reaction 

2AB^A 2 + B 2 , (3) 

the state A = corresponds to a molecule AB, whereas the state A = 1 corresponds to A 2 
+ B 2 . Representing the progress of the chemical reaction by a real- value parameter seems 
questionable when a mesoparticle stands for a single molecule. The progress variable should 
therefore rather be seen as some dissociation probability, or progress along some free energy 
profile. 

Since the model reaction (J3} involves two species on each side, we postulate for example 
the reversible evolution: 

= Y, u '-(ri i )[K 1 (T ij )(l-Xi)(l-X j )+K2{T ij )\ i X j \ (4) 

when < A; < 1, and ^ = otherwise (in order to ensure that the progress variables 
remains in [0, 1]). The function uj r in @ is a weight function (with support in [0,r rcac ]), and 
the mean temperature Tij is defined as in $2§. We term this reaction 'reversible' since Ai 
can either increase or decrease. The motivation for the postulated kinetics Q is the physical 
picture of a second-order reaction, where two molecules have to interact for the dissociation 
process to occur. Of course, many other kinetics are also possible, such as the first-order 
irreversible evolution 

where (h) = Y2j u> r (rij)hi denotes a local spatial average. The choice of the reaction kinetics 
is really a modelling choice, depending on the physical context. 

The reaction constants K\, K% are assumed to depend only on internal temperatures of 
particles. For example, a possible form in the Arrhenius spirit is: 

Ki(T) = Zx e~ El/kBT , K 2 (T) = Z 2 e - E2/kBT . (5) 

In this expression, Ex and £" 2 are the activation energies of, respectively, the forward and 
backward reactions. In view of our choice Q of reaction kinetics, the total increment of the 
progress variable is therefore the sum of all elementary pair increments, which is very much 
in the DPD spirit . The parameters Zx , Zi , Ex , E% in the above equation can be obtained from 
all-atom decomposition simulations or directly from experiments. In particular, Zx and Ex can 
be extracted in the Arrhenius framework from the analysis of decomposition rates at different 
initial temperatures. The constants Z 2 and E2 can be evaluated using thermochemical data 
(formation enthalpies). 

For very exothermic reactions, E2 3> Ex, and both energies are large since the activation 
energy is usually large for energetic materials. The increment of a given reaction rate is 
therefore non-negligible only if the material is sufficiently heated. In practice, this can be 
achieved when a strong shock is initiated in the system. If this shock is not strong enough, 
chemical reactions do not occur fast enough, and since the energy release is not sufficient, 
the shock wave is weakened until it transforms into a sonic wave. On the contrary, if the 
shock wave is strong enough, chemical reactions happen close enough to the detonation front, 
and the energy released sustains the shock wave such that it transforms into a stationary 
detonation wave. The progress of the reaction also modifies the mechanical properties of the 
material. In particular, reaction products usually have a larger specific volume than reactants 
(at fixed thermodynamic conditions). Therefore, some expansion is expected. The changes in 
the nature of the molecules are taken into account by introducing two additional parameters 
k a ,kE and using some mixing rule such as Berthelot's rule. When the interaction potential 
is of Lennard- Jones form, the interaction between the mesoparticles i and j is then given by 

(6) 
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with 

Eij = E^(l + k E \i){\ + k E \j), 

and 

A A; + Aj \ 

ai-j = a I 1 + k a I . 

When the reaction is complete, the material initially described by a Lennard- Jones potential 
of parameters a,E is then described by a Lennard- Jones of parameters a' = a(l + k a ) and 
E' = E(l + k E ). 



4 Treatment of the exothermicity 

We denote by A-Eexthm the exothermicity of the reaction (J3}. It is expected that A_E cxt hm = 
Ei — Ei. We assume that the energy is liberated progressively during the reaction, in a 
manner that the total energy of the system (chemical, mechanical, internal) is preserved: 



dHt, 



i=l 1 l<i<j<AT 



0. 



In order to propose a dynamics satisfying this condition, we have to make an additional 
assumption about the evolution of the system. Neglecting diffusive processes, we require 
that, during the elementary step corresponding to exothermicity, the total energy of a given 
mesoparticle does not changj]]: 



, Ai , \j ) 



+ d 



2mi 



AE . 



Xi = o. 



(7) 



We then consider evolutions of momenta and internal energies balancing the variations in the 
total energy due to the variations of A (exothermicity, changes in potential energies). This is 
analogous to the fact that the variations of kinetic energy in ([TJ) are compensated by varia- 
tions of internal energies. It is expected that the chemical energy liberated by the reaction is 
converted into internal energy and kinetic energy. It therefore remains to precise the quanti- 
tative distribution of the energy. This is done using a predetermined ratio < c < 1. This 
ratio could be measured in gas phase decomposition experiments, using spectroscopy mea- 
surements (possibly numerical simulations as well, by computing the temperature of internal 
degrees of freedom after the reaction). Alternatively, it is possible to postulate that c should 
be comparable to the ratio of the number of external degrees of freedom for the products of 
the dissociation of one complex molecule (denoted by N c ), divided by the total number of 
degrees of freedom N t in a complex molecule (assuming some equipartition of the liberated 
chemical energy). For exothermic chemical decompositions, it is expected that a single com- 
plex molecule breaks into several smaller molecules (still aggregated in a single mesoparticle), 
so that the above ratio N c /N t should be close to 0.5. Let us however emphasize at this point 
that numerical profiles obtained here are robust enough with respect to the choice of c. 

For the internal energies, the exothermicity of the reaction is modelled as dti = d£i with 



d£i 



-c d 



AE C 



For the momenta, we consider process V% with dpi = dVi such that 



-(l-c) d 



i^V^A^Aj 



a dXi 



We explain in the next section how this is done in practice (see Eq. (TlOlO . 

Let us emphasize at this point that there are many other possible ways to treat the 
exothermicity. For instance, it would be possible to consider instantaneous reactions (jump 



x Of course, during the elementary step corresponding to the dynamics {TJ, the total energy changes. 
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processes for which A changes from to 1) occuring at random times, the probability of 
reaction depending on the progress variable. However, it is not clear whether such a dynamics 
is reversible, since the reverse reaction requires particles to have large kinetic and internal 
energies. In comparison, the process described here is progressive and therefore, much more 
reversible. 

Finally, we propose the following dynamics to describe reactive shock waves: 

dqt = — dt, 

rrii 

dpi = 22 -V 9i V(nj, Xi, Xj) dt - fijX (rij)vij dt + ax(rij)dWij + dVi, (8) 

dti = ^ X 2 (fij) (jijVij - ^- [~ + ~)) dt -° X(nj)vij ■ dWij + dSi, 
dX t = ^wr(rij) [K^Xl - Xi){l - Xj) + K 2 (Tij)XiXj] dt, 

where dVi, d£i are such that (JT} holds, i.e. the total energy is conserved. The fluctuation- 
dissipation relation relating 7y and a is still ([2]). Notice also that the inert dynamics (HJ of [9] 
is recovered when Z\ = Zi = 0, starting from Xi = for all i. 



5 Numerical implementation 



dXi = 




dpi = 


dVi, 


dei = 


d£i. 



+ ^>,(rr 3 +1 ) [Kx{f^){l - A?)(l - A?)tf 2 (7$)A?A?] At. 



The numerical integration of (O is done using a decomposition of the dynamics into ele- 
mentary stochastic differential equations. We denote by <^>i ner t the flow associated with the 
dynamics (!]), and by </>*eao the flow associated with the remaining part of the dynamics JU): 

VI < i < AT, 

{ dX z = E,^^(r IJ )[A'i(r i:J )(l-AO(l-A J ) + A' 2 (T i3 )A l A J ] dt, 

(9) 

A one-step integrator for a time-step At is constructed as 

/ n+1 n + 1 n + 1 \ n+l\ i At iAt / n n n \ n \ 

(<? ,P ,e ,X ) = $ reac o ® inert (g ,p ,e ,A ). 

A possible numerical flow <3?^ert is given in [9]. 

Let us now construct a numerical flow $^i c approximating the flow 0reac Denoting 
(q n+1 ,p n , e™, A") = $^e rt (q n ,p n , e", A"), we first integrate the evolution equation on the 
progress variables Xi using a first-order explicit integration: 

A™ +1 = Ar 

We then set A™ +1 = min(max(0, A™ +1 ), 1) in order to ensure that the progress variable remains 
between and 1. Once all progress variables are updated, the variation 5E™ in the total energy 
of particle i due to the variations of {Xj } is computed as 

8Ef = (Ar +1 - X?)AE cxthm + i J2 {V(rZ +1 ,X?+\X] +1 ) - V(r^\X?, A™)) . 

The conservation of total energy is then ensured through variations of internal and kinetic 
energies. The internal energies are updated as = e™ + c5E™. The update of p" +1 is done 
by adding to p™ a vector with random direction, so that the final momentum is such that the 
kinetic energy is correct. More precisely, when the dimension of the physical space is d = 2 
for example, an angle 9f is chosen at random in the interval [0, 2ty], the angles (#™)j,n being 
independent and identically distributed random variables. The new momentum is then 
constructed such that 

n+1 n n/ r\n • /in\ 

Pi =Pi +a (cos0i ,sin0i ), 

2rrn 2rm y ' 1 y ' 

Solving this equation in a n gives the desired result. 
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6 Numerical results 



We present in this section numerical results obtained for the dynamics JHJ of a two-dimensional 
fluid. A shock is initiated using a piston of velocity u p during a time t p . The initial conditions 
for the positions qi, momenta Pi and internal energies are sampled as proposed in [§]■ 

We consider the following parameters, inspired by the nitromethane example (in which 
case the molecule CH3NO2 would replaced by a mesoparticle in a space of 2 dimensions). 
The parameters can be classified in four main categories, depending on whether they describe 
the mechanical properties of the material, characterize the inert dynamics and the chemical 
kinetics, or are related to the exothermicity. We consider here a system with 

• (Material parameters) a molar mass m = 80 g/mol, described by a Lennard- Jones 
potential of parameter £"lj = 3 x 1CP 21 J (melting temperature around 220 K) and 
a — 5 A, with a cut-off radius r cut = 15 A for the computation of forces. The changes 
in the parameters of the Lennard- Jones material during the reaction follow J6}, using 
A:_b = and k a = 0.2 (pure expansion). 

• (Parameters of the inert dynamics) The microscopic state law is e = C V T with C v = 
10 ks (i.e., 20 degrees of freedom are not represented). The friction is 7 = 10~ 15 kg/s, 
and the dissipation weighting function x(f) = (1 — r/r c ), with r c = r cut . 

• (Chemical kinetics) For the chemical reaction ([4]), reaction constants are computed 
using (|5j) with Z\ — Z 2 — 10 17 s -1 , Ei/fc B = 15000 K, the exothermicity being 
A-Eexthm = 6.25 eV. The reaction weighting function w(r) = x( r )i 

• (Exothermicity) we choose c = 0.5. 

The initial density of the system is p = 1.06 g/cm -3 , and the initial temperature T — 300 K. 
The time-step used is At = 2 x 10 -15 s. Figure [l] presents velocity profiles averaged in thin 
slices of the material in the direction of the shock, for a compression time t p = 2 ps at a velocity 
Up = 5000 m/s. We tested the independence of velocity profiles in the reaction zone with 
respect to different initial loadings (t p ,u p ) = (1 ps, 6000 m/s), (t p ,u p ) — (2 ps, 6000 m/s), 
(tp,u p ) = (3 ps, 6000 m/s) and (t p ,u p ) = (3 ps, 5000 m/s). 

The velocity of the shock front is constant, and approximately equal to u s = 1820 m/s. 
Notice that the wave can be divided into three regions: the upstream region is unperturbed; 
the region around the shock front where chemical reactions happen is of constant width 
(approximately 300-400 A, which is consistent with all-atoms studies, see for instance |16|): 
the downstream region is an autosimilar rarefaction wave. This profile is therefore reminiscent 
of ZND profiles [T] encountered in hydrodynamic simulations of detonation waves. 

Figure [2] presents the evolution of internal and kinetic temperatures averaged in a slice of 
material in the direction of the shock as a function of time (Left), as well as the evolution of 
the average progress variables (Right). In particular, the reaction does not start immedialely 
at the shock front: the ignition asks first for a sufficient heating of the material (through 
an increasing internal energy), since the reaction constant are too low at temperatures lower 
than a few thousands Kelvins with the values chosen here. 

7 Conclusion and perspectives 

In conclusion, we extended the mesoscopic inert model for shock waves of [9] to the reactive 
case. The key idea is to introduce an additional variable describing the progress of the reaction, 
and to postulate an activated exothermic chemical reaction. Provided the initial loading is 
strong enough, the energy released by chemical reactions sustains the shock wave, whose front 
then travels unperturbed in the material, followed by a rarefaction wave. 

Now that while the qualitative behavior of the model is granted, more quantitative studies 
must be pursued, where a careful comparison between all-atom numerical results and the 
profiles predicted by the reduced model proposed here are to be compared. Such studies were 
already performed in the inert case [8j [9] . 

A promising track for a real multiscale modelling of detonation waves would now be 
to further reduce the model proposed here, in order to obtain a description of matter at 
a scale truly intermediate between SPH and all-atom models. In such a model, random 
fluctuations would still be present, though of much lower magnitude; on the other hand, local 



6 



3000 




-0.4 -0.3 -0.2 -0.1 0.0 0.1 

Distance to the front (\im) 

Figure 1: (color online) Velocity profiles in the material at different times, with shock fronts 
aligned (Black: t = 1.6 x KT 10 s; Blue: t = 2.0 x 1(T 10 s; Red: t = 2.4 x 1CT 10 s). 
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Figure 2: (color online) Left : variations of internal (red) and kinetic (black) temperatures in the 
direction of the shock, as a function of time in a slice of material. Right: evolution of the progress 
variable averaged in a slice of material as a function of time (blue). For comparison, a reseated 
internal temperature profile is also presented (black). 
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thermodynamic fields (such as pressure) should also be introduced. Such an intermediate 
model is necessary to describe mesoscopic effects having an impact on the macroscopic features 
of the detonation, such as hot spot ignition. 
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